<!DOCTYPE html PUBLIC "-//W3C//DTD HTML 4.01 Transitional//EN">
<html><head><script type="text/javascript" src="https://cdn.jsdelivr.net/gh/opencobra/cobratoolbox@ffa0229fc0c01c9236bb7e961f65712443277719/latest/_static/js/iframeResizer.contentWindow.min.js"></script><meta http-equiv="Content-Type" content="text/html; charset=utf-8"><meta http-equiv="X-UA-Compatible" content="IE=edge,IE=9,chrome=1"><meta name="generator" content="MATLAB 2021a"><title>Sensitivity of a flux balance analysis solution with respect to input data</title><style type="text/css">.rtcContent { padding: 30px; } .S0 { margin: 3px 10px 5px 4px; padding: 0px; line-height: 28.8px; min-height: 0px; white-space: normal; color: rgb(213, 80, 0); font-family: Helvetica, Arial, sans-serif; font-style: normal; font-size: 24px; font-weight: normal; text-align: left;  }
.S1 { margin: 20px 10px 5px 4px; padding: 0px; line-height: 20px; min-height: 0px; white-space: normal; color: rgb(60, 60, 60); font-family: Helvetica, Arial, sans-serif; font-style: normal; font-size: 20px; font-weight: bold; text-align: left;  }
.S2 { margin: 2px 10px 9px 4px; padding: 0px; line-height: 21px; min-height: 0px; white-space: normal; color: rgb(0, 0, 0); font-family: Helvetica, Arial, sans-serif; font-style: normal; font-size: 14px; font-weight: normal; text-align: left;  }
.CodeBlock { background-color: #F7F7F7; margin: 10px 0 10px 0;}
.S3 { border-left: 1px solid rgb(233, 233, 233); border-right: 1px solid rgb(233, 233, 233); border-top: 1px solid rgb(233, 233, 233); border-bottom: 0px none rgb(0, 0, 0); border-radius: 4px 4px 0px 0px; padding: 6px 45px 0px 13px; line-height: 17.234px; min-height: 18px; white-space: nowrap; color: rgb(0, 0, 0); font-family: Menlo, Monaco, Consolas, "Courier New", monospace; font-size: 14px;  }
.S4 { border-left: 1px solid rgb(233, 233, 233); border-right: 1px solid rgb(233, 233, 233); border-top: 0px none rgb(0, 0, 0); border-bottom: 0px none rgb(0, 0, 0); border-radius: 0px; padding: 0px 45px 0px 13px; line-height: 17.234px; min-height: 18px; white-space: nowrap; color: rgb(0, 0, 0); font-family: Menlo, Monaco, Consolas, "Courier New", monospace; font-size: 14px;  }
.S5 { border-left: 1px solid rgb(233, 233, 233); border-right: 1px solid rgb(233, 233, 233); border-top: 0px none rgb(0, 0, 0); border-bottom: 1px solid rgb(233, 233, 233); border-radius: 0px 0px 4px 4px; padding: 0px 45px 4px 13px; line-height: 17.234px; min-height: 18px; white-space: nowrap; color: rgb(0, 0, 0); font-family: Menlo, Monaco, Consolas, "Courier New", monospace; font-size: 14px;  }
.S6 { margin: 10px 10px 9px 4px; padding: 0px; line-height: 21px; min-height: 0px; white-space: normal; color: rgb(0, 0, 0); font-family: Helvetica, Arial, sans-serif; font-style: normal; font-size: 14px; font-weight: normal; text-align: center;  }
.S7 { margin: 10px 0px 20px; padding-left: 0px; font-family: Helvetica, Arial, sans-serif; font-size: 14px;  }
.S8 { margin-left: 56px; line-height: 21px; min-height: 0px; text-align: left; white-space: normal;  }
.S9 { margin: 10px 10px 5px 4px; padding: 0px; line-height: 18px; min-height: 0px; white-space: normal; color: rgb(60, 60, 60); font-family: Helvetica, Arial, sans-serif; font-style: normal; font-size: 15px; font-weight: bold; text-align: left;  }
.S10 { border-left: 1px solid rgb(233, 233, 233); border-right: 1px solid rgb(233, 233, 233); border-top: 0px none rgb(0, 0, 0); border-bottom: 1px solid rgb(233, 233, 233); border-radius: 0px; padding: 0px 45px 4px 13px; line-height: 17.234px; min-height: 18px; white-space: nowrap; color: rgb(0, 0, 0); font-family: Menlo, Monaco, Consolas, "Courier New", monospace; font-size: 14px;  }
.S11 { color: rgb(64, 64, 64); padding: 10px 0px 6px 17px; background: rgb(255, 255, 255) none repeat scroll 0% 0% / auto padding-box border-box; font-family: Menlo, Monaco, Consolas, "Courier New", monospace; font-size: 14px; overflow-x: hidden; line-height: 17.234px;  }
/* Styling that is common to warnings and errors is in diagnosticOutput.css */.embeddedOutputsErrorElement {    min-height: 18px;    max-height: 250px;    overflow: auto;}
.embeddedOutputsErrorElement.inlineElement {}
.embeddedOutputsErrorElement.rightPaneElement {}
/* Styling that is common to warnings and errors is in diagnosticOutput.css */.embeddedOutputsWarningElement{    min-height: 18px;    max-height: 250px;    overflow: auto;}
.embeddedOutputsWarningElement.inlineElement {}
.embeddedOutputsWarningElement.rightPaneElement {}
/* Copyright 2015-2019 The MathWorks, Inc. *//* In this file, styles are not scoped to rtcContainer since they could be in the Dojo Tooltip */.diagnosticMessage-wrapper {    font-family: Menlo, Monaco, Consolas, "Courier New", monospace;    font-size: 12px;}
.diagnosticMessage-wrapper.diagnosticMessage-warningType {    color: rgb(255,100,0);}
.diagnosticMessage-wrapper.diagnosticMessage-warningType a {    color: rgb(255,100,0);    text-decoration: underline;}
.diagnosticMessage-wrapper.diagnosticMessage-errorType {    color: rgb(230,0,0);}
.diagnosticMessage-wrapper.diagnosticMessage-errorType a {    color: rgb(230,0,0);    text-decoration: underline;}
.diagnosticMessage-wrapper .diagnosticMessage-messagePart,.diagnosticMessage-wrapper .diagnosticMessage-causePart {    white-space: normal;}
.diagnosticMessage-wrapper .diagnosticMessage-stackPart {    white-space: normal;}
.embeddedOutputsTextElement,.embeddedOutputsVariableStringElement {    white-space: normal;    word-wrap:  initial;    min-height: 18px;    max-height: 250px;    overflow: auto;}
.textElement,.rtcDataTipElement .textElement {    padding-top: 3px;}
.embeddedOutputsTextElement.inlineElement,.embeddedOutputsVariableStringElement.inlineElement {}
.inlineElement .textElement {}
.embeddedOutputsTextElement.rightPaneElement,.embeddedOutputsVariableStringElement.rightPaneElement {    min-height: 16px;}
.rightPaneElement .textElement {    padding-top: 2px;    padding-left: 9px;}
.S12 { border-left: 1px solid rgb(233, 233, 233); border-right: 1px solid rgb(233, 233, 233); border-top: 1px solid rgb(233, 233, 233); border-bottom: 0px none rgb(0, 0, 0); border-radius: 0px; padding: 6px 45px 0px 13px; line-height: 17.234px; min-height: 18px; white-space: nowrap; color: rgb(0, 0, 0); font-family: Menlo, Monaco, Consolas, "Courier New", monospace; font-size: 14px;  }
.S13 { margin: 10px 10px 9px 4px; padding: 0px; line-height: 21px; min-height: 0px; white-space: normal; color: rgb(0, 0, 0); font-family: Helvetica, Arial, sans-serif; font-style: normal; font-size: 14px; font-weight: normal; text-align: left;  }
.S14 { border-left: 1px solid rgb(233, 233, 233); border-right: 1px solid rgb(233, 233, 233); border-top: 1px solid rgb(233, 233, 233); border-bottom: 1px solid rgb(233, 233, 233); border-radius: 4px 4px 0px 0px; padding: 6px 45px 4px 13px; line-height: 17.234px; min-height: 18px; white-space: nowrap; color: rgb(0, 0, 0); font-family: Menlo, Monaco, Consolas, "Courier New", monospace; font-size: 14px;  }
.S15 { border-left: 1px solid rgb(233, 233, 233); border-right: 1px solid rgb(233, 233, 233); border-top: 1px solid rgb(233, 233, 233); border-bottom: 1px solid rgb(233, 233, 233); border-radius: 0px; padding: 6px 45px 4px 13px; line-height: 17.234px; min-height: 18px; white-space: nowrap; color: rgb(0, 0, 0); font-family: Menlo, Monaco, Consolas, "Courier New", monospace; font-size: 14px;  }</style></head><body><div class = rtcContent><h1  class = 'S0'><span>Sensitivity of a flux balance analysis solution with respect to input data</span></h1><h2  class = 'S1'><span>Author(s): Ronan M.T. Fleming, Leiden University</span></h2><h2  class = 'S1'><span>Reviewer(s): </span></h2><h2  class = 'S1'><span>INTRODUCTION</span></h2><div  class = 'S2'><span>Consider an FBA problem</span></div><div  class = 'S2'><span texencoding="
\begin{array}{ll}
\textrm{max} &amp; c^{T}v\\
\text{s.t.} &amp; Sv=b\\
 &amp; l\leq v\leq u
\end{array}
\end{equation}" style="vertical-align:-25px"><img src="" width="96.5" height="61" /></span></div><div  class = 'S2'><span>The local sensitivity of the optimal objective value </span><span texencoding="\mathcal{L}^{\star} = c^{T}v^{\star}" style="vertical-align:-5px"><img src="" width="69.5" height="19" /></span><span>  with respect to a changes in the input data </span><span mathmlencoding="&lt;math xmlns=&quot;http://www.w3.org/1998/Math/MathML&quot; display=&quot;inline&quot;&gt;&lt;mrow&gt;&lt;mrow&gt;&lt;mo&gt;{&lt;/mo&gt;&lt;mrow&gt;&lt;mi mathvariant=&quot;italic&quot;&gt;b&lt;/mi&gt;&lt;mo&gt;,&lt;/mo&gt;&lt;mi mathvariant=&quot;italic&quot;&gt;l&lt;/mi&gt;&lt;mo&gt;,&lt;/mo&gt;&lt;mi mathvariant=&quot;italic&quot;&gt;u&lt;/mi&gt;&lt;/mrow&gt;&lt;mo&gt;}&lt;/mo&gt;&lt;/mrow&gt;&lt;/mrow&gt;&lt;/math&gt;" style="vertical-align:-5px"><img src="" width="50" height="19" /></span><span> is given by</span></div><div  class = 'S2'><span> </span><span texencoding="\frac{\partial\mathcal{L}^{\star}}{\partial b} = y^{\star} \\
\frac{\partial\mathcal{\mathcal{L}^{\star}}}{\partial l} = -w_{l}^{\star} \\
\frac{\partial\mathcal{\mathcal{L}^{\star}}}{\partial u} = w_{u}^{\star}
" style="vertical-align:-53px"><img src="" width="77" height="117" /></span></div><div  class = 'S2'><span>where </span><span texencoding="y^{\star}" style="vertical-align:-5px"><img src="" width="17.5" height="19" /></span><span> is a vector of shadow prices and </span><span texencoding="w = w_{l}-w_{u}" style="vertical-align:-6px"><img src="" width="76.5" height="20" /></span><span> is a vector of reduced costs. That is, a shadow price is the partial derivative of the optimal value of the objective function with respect to </span><span mathmlencoding="&lt;math xmlns=&quot;http://www.w3.org/1998/Math/MathML&quot; display=&quot;inline&quot;&gt;&lt;mrow&gt;&lt;msub&gt;&lt;mrow&gt;&lt;mi mathvariant=&quot;italic&quot;&gt;b&lt;/mi&gt;&lt;/mrow&gt;&lt;mrow&gt;&lt;mi mathvariant=&quot;italic&quot;&gt;i&lt;/mi&gt;&lt;/mrow&gt;&lt;/msub&gt;&lt;/mrow&gt;&lt;/math&gt;" style="vertical-align:-6px"><img src="" width="13.5" height="20" /></span><span>.  It indicates how much net production, or net consumption, of each metabolite increases (positive), or decreases (negative), the optimal value of the objective. The reduced costs, </span><span mathmlencoding="&lt;math xmlns=&quot;http://www.w3.org/1998/Math/MathML&quot; display=&quot;inline&quot;&gt;&lt;mrow&gt;&lt;mo&gt;-&lt;/mo&gt;&lt;msub&gt;&lt;mrow&gt;&lt;mi mathvariant=&quot;italic&quot;&gt;w&lt;/mi&gt;&lt;/mrow&gt;&lt;mrow&gt;&lt;mi mathvariant=&quot;italic&quot;&gt;l&lt;/mi&gt;&lt;/mrow&gt;&lt;/msub&gt;&lt;/mrow&gt;&lt;/math&gt;" style="vertical-align:-6px"><img src="" width="26" height="20" /></span><span> and </span><span mathmlencoding="&lt;math xmlns=&quot;http://www.w3.org/1998/Math/MathML&quot; display=&quot;inline&quot;&gt;&lt;mrow&gt;&lt;msub&gt;&lt;mrow&gt;&lt;mi mathvariant=&quot;italic&quot;&gt;w&lt;/mi&gt;&lt;/mrow&gt;&lt;mrow&gt;&lt;mi mathvariant=&quot;italic&quot;&gt;u&lt;/mi&gt;&lt;mtext&gt; &lt;/mtext&gt;&lt;/mrow&gt;&lt;/msub&gt;&lt;/mrow&gt;&lt;/math&gt;" style="vertical-align:-6px"><img src="" width="20" height="20" /></span><span>are the partial derivative of the optimal value of the objective function with respect to the lower and upper bounds on a reaction, respectively. They indicate how much relaxation, or tightening, of each bound increases, or decreases, the optimal objective, respectively. In the COBRA Toolbox, shadow prices and reduced costs are calculated by </span><span style=' font-family: monospace;'>optimizeCbModel</span><span>. When using the function</span></div><div  class = 'S2'><span style=' font-family: monospace;'>FBAsolution = optimizeCbModel(model,'max');</span></div><div  class = 'S2'><span>the shadow prices and reduced costs are given by </span><span style=' font-family: monospace;'>FBAsolution.y </span><span>and</span><span style=' font-family: monospace;'> FBAsolution.w</span><span>, respectively. </span></div><div  class = 'S2'><span>For a more complete theoretical description, see: cobratoolbox/tutorials/intro_sensitivityAnalysis.pdf</span></div><h2  class = 'S1'><span>MATERIALS - EQUIPMENT SETUP</span></h2><div  class = 'S2'><span>Please ensure that all the required dependencies (e.g. , </span><span style=' font-family: monospace;'>git</span><span> and </span><span style=' font-family: monospace;'>curl</span><span>) of The COBRA Toolbox have been properly installed by following the installation guide </span><a href = "https://opencobra.github.io/cobratoolbox/stable/installation.html"><span>here</span></a><span>. Please ensure that the COBRA Toolbox has been initialised (tutorial_initialize.mlx) and verify that the pre-packaged LP and QP solvers are functional (tutorial_verify.mlx).</span></div><h2  class = 'S1'><span>PROCEDURE</span></h2><h2  class = 'S1'><span>Load E. coli core model</span></h2><div  class = 'S2'><span>The most direct way to load a model into The COBRA Toolbox is to use the </span><span style=' font-family: monospace;'>readCbModel</span><span> function. For example, to load a model from a MAT-file, you can simply use the filename (with or without file extension). </span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S3'><span style="white-space: normal"><span >fileName = </span><span style="color: rgb(170, 4, 249);">'ecoli_core_model.mat'</span><span >;</span></span></div></div><div class="inlineWrapper"><div  class = 'S4'><span style="white-space: normal"><span style="color: rgb(14, 0, 255);">if </span><span >~exist(</span><span style="color: rgb(170, 4, 249);">'modelOri'</span><span >,</span><span style="color: rgb(170, 4, 249);">'var'</span><span >)</span></span></div></div><div class="inlineWrapper"><div  class = 'S4'><span style="white-space: normal"><span >modelOri = readCbModel(fileName);</span></span></div></div><div class="inlineWrapper"><div  class = 'S4'><span style="white-space: normal"><span style="color: rgb(14, 0, 255);">end</span></span></div></div><div class="inlineWrapper"><div  class = 'S4'><span style="white-space: normal"><span style="color: rgb(2, 128, 9);">%backward compatibility with primer requires relaxation of upper bound on</span></span></div></div><div class="inlineWrapper"><div  class = 'S4'><span style="white-space: normal"><span style="color: rgb(2, 128, 9);">%ATPM</span></span></div></div><div class="inlineWrapper"><div  class = 'S4'><span style="white-space: normal"><span >modelOri = changeRxnBounds(modelOri,</span><span style="color: rgb(170, 4, 249);">'ATPM'</span><span >,1000,</span><span style="color: rgb(170, 4, 249);">'u'</span><span >);</span></span></div></div><div class="inlineWrapper"><div  class = 'S4'><span style="white-space: normal"><span >model = modelOri;</span></span></div></div><div class="inlineWrapper"><div  class = 'S4'><span style="white-space: normal"><span style="color: rgb(2, 128, 9);">%setp the matlab e.coli metabolic map parameters</span></span></div></div><div class="inlineWrapper"><div  class = 'S4'><span style="white-space: normal"><span >outputFormatOK = changeCbMapOutput(</span><span style="color: rgb(170, 4, 249);">'matlab'</span><span >);</span></span></div></div><div class="inlineWrapper"><div  class = 'S4'><span style="white-space: normal"><span >map=readCbMap(</span><span style="color: rgb(170, 4, 249);">'ecoli_core_map'</span><span >);</span></span></div></div><div class="inlineWrapper"><div  class = 'S4'><span style="white-space: normal"><span >options.zeroFluxWidth = 0.1;</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >options.rxnDirMultiplier = 10;</span></span></div></div></div><div  class = 'S6'><img class = "imageNode" src = "" alt = "" style = "vertical-align: baseline"></img></div><div  class = 'S2'><span>The meaning of each field in a standard model is defined in the </span><a href = "https://github.com/opencobra/cobratoolbox/blob/master/docs/source/notes/COBRAModelFields.md"><span>standard COBRA model field definition</span></a><span>.</span></div><div  class = 'S2'><span>In general, the following fields should always be present: </span></div><ul  class = 'S7'><li  class = 'S8'><span style=' font-weight: bold;'>S</span><span>, the stoichiometric matrix</span></li><li  class = 'S8'><span style=' font-weight: bold;'>mets</span><span>, the identifiers of the metabolites</span></li><li  class = 'S8'><span style=' font-weight: bold;'>b</span><span>, Accumulation (positive) or depletion (negative) of the corresponding metabolites. 0 Indicates no concentration change.</span></li><li  class = 'S8'><span style=' font-weight: bold;'>csense</span><span>, indicator whether the b vector is a lower bound ('G'), upper bound ('L'), or hard constraint 'E' for the metabolites.</span></li><li  class = 'S8'><span style=' font-weight: bold;'>rxns</span><span>, the identifiers of the reactions</span></li><li  class = 'S8'><span style=' font-weight: bold;'>lb</span><span>, the lower bounds of the reactions</span></li><li  class = 'S8'><span style=' font-weight: bold;'>ub</span><span>, the upper bounds of the reactions</span></li><li  class = 'S8'><span style=' font-weight: bold;'>c</span><span>, the linear objective</span></li><li  class = 'S8'><span style=' font-weight: bold;'>genes</span><span>, the list of genes in your model </span></li><li  class = 'S8'><span style=' font-weight: bold;'>rules</span><span>, the Gene-protein-reaction rules in a computer readable format present in your model.</span></li><li  class = 'S8'><span style=' font-weight: bold;'>osenseStr</span><span>, the objective sense either </span><span style=' font-family: monospace;'>'max'</span><span> for maximisation or </span><span style=' font-family: monospace;'>'min'</span><span> for minimisation</span></li></ul><h2  class = 'S1'><span>Sensitivity Analysis</span></h2><h4  class = 'S9'><span>In the E. coli core model, when maximising ATP production, what is the shadow price of cytosolic protons? </span></h4><h4  class = 'S9'><span>Hint: </span><span style=' font-family: monospace;'>FBAsolution.y</span></h4><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S3'><span style="white-space: normal"><span >model = modelOri;</span></span></div></div><div class="inlineWrapper"><div  class = 'S4'><span style="white-space: normal"><span >model = changeRxnBounds(model,</span><span style="color: rgb(170, 4, 249);">'EX_glc(e)'</span><span >,-1,</span><span style="color: rgb(170, 4, 249);">'l'</span><span >);</span></span></div></div><div class="inlineWrapper"><div  class = 'S4'><span style="white-space: normal"><span >model = changeRxnBounds(model,</span><span style="color: rgb(170, 4, 249);">'EX_o2(e)'</span><span >,-1000,</span><span style="color: rgb(170, 4, 249);">'l'</span><span >);</span></span></div></div><div class="inlineWrapper"><div  class = 'S4'><span style="white-space: normal"><span >model = changeRxnBounds(model,</span><span style="color: rgb(170, 4, 249);">'ATPM'</span><span >,0,</span><span style="color: rgb(170, 4, 249);">'l'</span><span >);</span></span></div></div><div class="inlineWrapper"><div  class = 'S4'><span style="white-space: normal"><span >model = changeObjective(model,</span><span style="color: rgb(170, 4, 249);">'ATPM'</span><span >);</span></span></div></div><div class="inlineWrapper outputs"><div  class = 'S10'><span style="white-space: normal"><span >printConstraints(model,-1000,1000)</span></span></div><div  class = 'S11'><div class="inlineElement eoOutputWrapper embeddedOutputsTextElement" uid="2FA01D15" data-testid="output_0" data-width="428" data-height="45" data-hashorizontaloverflow="false" style="width: 458px; max-height: 261px; white-space: normal; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"><div class="textElement" style="white-space: normal; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;">MinConstraints:
EX_glc(e)	-1
maxConstraints:</div></div></div></div><div class="inlineWrapper"><div  class = 'S12'></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >FBAsolution_maxATP = optimizeCbModel(model,</span><span style="color: rgb(170, 4, 249);">'max'</span><span >);</span></span></div></div></div><div  class = 'S13'><span>Check the optimal value of the objective</span></div><div class="CodeBlock"><div class="inlineWrapper outputs"><div  class = 'S14'><span style="white-space: normal"><span >FBAsolution_maxATP.f</span></span></div><div  class = 'S11'><div class='variableElement' style='font-family: Menlo, Monaco, Consolas, "Courier New", monospace; font-size: 12px; '>ans = 17.5000</div></div></div></div><div  class = 'S13'><span>The shadow price of cytosolic protons (h[c]) is -0.25. </span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S3'><span style="white-space: normal"><span >ind=strcmp(model.mets,</span><span style="color: rgb(170, 4, 249);">'h[c]'</span><span >);</span></span></div></div><div class="inlineWrapper outputs"><div  class = 'S10'><span style="white-space: normal"><span >FBAsolution_maxATP.y(ind)</span></span></div><div  class = 'S11'><div class='variableElement' style='font-family: Menlo, Monaco, Consolas, "Courier New", monospace; font-size: 12px; '>ans = -0.2500</div></div></div><div class="inlineWrapper outputs"><div  class = 'S15'><span style="white-space: normal"><span >printFluxVector(model, FBAsolution_maxATP.v, 1)</span></span></div><div  class = 'S11'><div class="inlineElement eoOutputWrapper embeddedOutputsTextElement" uid="0112DCAF" data-testid="output_3" data-width="428" data-height="451" data-hashorizontaloverflow="false" style="width: 458px; max-height: 462px; white-space: normal; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"><div class="textElement" style="white-space: normal; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;">ACONTa              	           2
ACONTb              	           2
AKGDH               	           2
ATPM                	        17.5
ATPS4r              	        13.5
CO2t                	          -6
CS                  	           2
CYTBD               	          12
ENO                 	           2
EX_co2(e)           	           6
EX_glc(e)           	          -1
EX_h2o(e)           	           6
EX_o2(e)            	          -6
FBA                 	           1
FUM                 	           2
GAPD                	           2
GLCpts              	           1
H2Ot                	          -6
ICDHyr              	           2
MDH                 	           2
NADH16              	          10
NADTRHD             	           2
O2t                 	           6
PDH                 	           2
PFK                 	           1
PGI                 	           1
PGK                 	          -2
PGM                 	          -2
PYK                 	           1
SUCDi               	           2
SUCOAS              	          -2
TPI                 	           1</div></div></div></div></div><h4  class = 'S9'><span>What is your biochemical interpretation of this change in objective in the current context?</span></h4><h4  class = 'S9'><span>Hint: printFluxVector, drawFlux</span></h4><div  class = 'S2'><span>This is a unique solution (see Example 3). </span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S3'><span style="white-space: normal"><span >dv = FBAsolution_maxATP_forceH.v-FBAsolution_maxATP.v;</span></span></div></div><div class="inlineWrapper"><div  class = 'S4'><span style="white-space: normal"><span >dv(abs(dv)&lt;1e-5)=0;</span></span></div></div><div class="inlineWrapper outputs"><div  class = 'S10'><span style="white-space: normal"><span >printFluxVector(model, dv, 1)</span></span></div><div  class = 'S11'><div class="inlineElement eoOutputWrapper embeddedOutputsTextElement" uid="16A81FB8" data-testid="output_4" data-width="428" data-height="45" data-hashorizontaloverflow="false" style="width: 458px; max-height: 261px; white-space: normal; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"><div class="textElement" style="white-space: normal; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;">ATPM                	           1
ATPS4r              	           1
EX_h(e)             	          -4</div></div></div></div></div><div  class = 'S13'><span>The flux map for optimal ATP production is shown below.  </span></div><div class="CodeBlock"><div class="inlineWrapper outputs"><div  class = 'S14'><span style="white-space: normal"><span >drawFlux(map, model, FBAsolution_maxATP.v, options);</span></span></div><div  class = 'S11'><div class="inlineElement eoOutputWrapper embeddedOutputsFigure" uid="E6C59849" data-testid="output_5" style="width: 458px;"><div class="figureElement"><img class="figureImage figureContainingNode" src="" style="width: 560px;"></div></div></div></div></div><div  class = 'S13'><span>ATP production is constrained by cytoplasmic proton balancing.  Cytoplasmic protons are produced by various metabolic reactions and also enter into the cell, from the extracellular compartment, via the ATP synthase reaction (ATPS4r).  At steady-state, an equal number of protons must be pumped out of the cytoplasm by the electron transport chain reactions or by excreting metabolites with symporters. Setting model.b(i) = 4, where i corresponds to cytoplasmic protons, h[c], removes 4 extra units of cytoplasmic protons from the system allowing 4 extra extracellular protons to enter the system that then enter the cell via the ATP synthase reaction, generating one extra unit of ATP. This increases the maximum rate of ATP synthesis by one unit, thereby increasing the ATP yield from glucose by 1 mol ATP/mol glucose. </span></div><h4  class = 'S9'><span>Perturb the model in such a way as to increase the optimal rate of ATP hydrolysis ('ATPM') by exactly one unit. How does this compare with the theoretical prediction?</span></h4><h4  class = 'S9'><span>Hint: change model.b</span></h4><div  class = 'S2'><span>Remove 4 units of cytoplasmic protons from the system, but changing model.b(i) to 4, where i corresponds to the index for cytoplasmic protons, and calculate the difference in the value of the optimal objective. The answer should be 1.</span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S3'><span style="white-space: normal"><span >ind=strcmp(model.mets,</span><span style="color: rgb(170, 4, 249);">'h[c]'</span><span >);</span></span></div></div><div class="inlineWrapper"><div  class = 'S4'><span style="white-space: normal"><span >model.b(ind) = 4;</span></span></div></div><div class="inlineWrapper"><div  class = 'S4'><span style="white-space: normal"><span >FBAsolution_maxATP_forceH = optimizeCbModel(model,</span><span style="color: rgb(170, 4, 249);">'max'</span><span >);</span></span></div></div><div class="inlineWrapper outputs"><div  class = 'S10'><span style="white-space: normal"><span >FBAsolution_maxATP_forceH.f - FBAsolution_maxATP.f</span></span></div><div  class = 'S11'><div class='variableElement' style='font-family: Menlo, Monaco, Consolas, "Courier New", monospace; font-size: 12px; '>ans = 1</div></div></div></div><h4  class = 'S9'><span>In the E. coli core model, when maximising ATP production, what is the reduced cost of glucose exchange? </span></h4><h4  class = 'S9'><span>Hint: FBAsolution.rcost</span></h4><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S3'><span style="white-space: normal"><span >rcost = FBAsolution_maxATP.rcost;</span></span></div></div><div class="inlineWrapper"><div  class = 'S4'><span style="white-space: normal"><span >rcost(abs(rcost)&lt;1e-4)=0;</span></span></div></div><div class="inlineWrapper"><div  class = 'S4'><span style="white-space: normal"><span >flux=FBAsolution_maxATP.v;</span></span></div></div><div class="inlineWrapper outputs"><div  class = 'S10'><span style="white-space: normal"><span >printFluxVector(model, [model.lb,flux,model.ub,rcost], 1)</span></span></div><div  class = 'S11'><div class="inlineElement eoOutputWrapper embeddedOutputsTextElement scrollableOutput" uid="E2D9AD08" data-testid="output_7" data-width="428" data-height="1333" data-hashorizontaloverflow="true" style="width: 458px; max-height: 261px; white-space: normal; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"><div class="textElement" style="white-space: normal; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;">ACALD               	       -1000           0        1000           0
ACALDt              	       -1000           0        1000           0
ACKr                	       -1000   1.202e-32        1000           0
ACONTa              	       -1000           2        1000           0
ACONTb              	       -1000           2        1000           0
ACt2r               	       -1000          -0        1000           0
ADK1                	       -1000           0        1000           0
AKGDH               	           0           2        1000           0
AKGt2r              	       -1000          -0        1000           0
ALCD2x              	       -1000           0        1000           0
ATPM                	           0        17.5        1000           0
ATPS4r              	       -1000        13.5        1000           0
Biomass_Ecoli_core_N(w/GAM)-Nmet2	           0           0        1000       188.3
CO2t                	       -1000          -6        1000           0
CS                  	           0           2        1000           0
CYTBD               	           0          12        1000           0
D-LACt2             	       -1000          -0        1000           0
ENO                 	       -1000           2        1000           0
ETOHt2r             	       -1000          -0        1000           0
EX_ac(e)            	           0           0        1000        4.25
EX_acald(e)         	           0           0        1000         6.5
EX_akg(e)           	           0           0        1000       11.75
EX_co2(e)           	       -1000           6        1000           0
EX_etoh(e)          	           0           0        1000         7.5
EX_for(e)           	           0          -0        1000           0
EX_fru(e)           	           0           0        1000        17.5
EX_fum(e)           	           0           0        1000        8.75
EX_glc(e)           	          -1          -1        1000        17.5
EX_gln-L(e)         	           0           0        1000       13.25
EX_glu-L(e)         	           0           0        1000          13
EX_h2o(e)           	       -1000           6        1000           0
EX_h(e)             	       -1000   1.449e-14        1000           0
EX_lac-D(e)         	           0           0        1000        7.75
EX_mal-L(e)         	           0           0        1000        8.75
EX_nh4(e)           	       -1000          -0        1000           0
EX_o2(e)            	       -1000          -6        1000           0
EX_pi(e)            	       -1000  -1.593e-16        1000           0
EX_pyr(e)           	           0           0        1000         6.5
EX_succ(e)          	           0           0        1000          10
FBA                 	       -1000           1        1000           0
FBP                 	           0           0        1000           1
FORt2               	           0           0        1000        0.25
FORti               	           0           0        1000           0
FRD7                	           0           0        1000           0
FRUpts2             	           0           0        1000           0
FUM                 	       -1000           2        1000           0
FUMt2_2             	           0           0        1000           0
G6PDH2r             	       -1000           0        1000           0
GAPD                	       -1000           2        1000           0
GLCpts              	           0           1        1000           0
GLNS                	           0           0        1000           0
GLNabc              	           0           0        1000           0
GLUDy               	       -1000           0        1000           0
GLUN                	           0           0        1000           1
GLUSy               	           0           0        1000           1
GLUt2r              	       -1000          -0        1000           0
GND                 	           0           0        1000      0.4167
H2Ot                	       -1000          -6        1000           0
ICDHyr              	       -1000           2        1000           0
ICL                 	           0           0        1000           0
LDH_D               	       -1000           0        1000           0
MALS                	           0           0        1000           0
MALt2_2             	           0           0        1000           0
MDH                 	       -1000           2        1000           0
ME1                 	           0           0        1000           1
ME2                 	           0           0        1000           1
NADH16              	           0          10        1000           0
NADTRHD             	           0           2        1000           0
NH4t                	       -1000           0        1000           0
O2t                 	       -1000           6        1000           0
PDH                 	           0           2        1000           0
PFK                 	           0           1        1000           0
PFL                 	           0           0        1000         1.5
PGI                 	       -1000           1        1000           0
PGK                 	       -1000          -2        1000           0
PGL                 	           0           0        1000           0
PGM                 	       -1000          -2        1000           0
PIt2r               	       -1000   1.593e-16        1000           0
PPC                 	           0  -3.403e-16        1000           0
PPCK                	           0           0        1000           1
PPS                 	           0           0        1000           1
PTAr                	       -1000  -1.202e-32        1000           0
PYK                 	           0           1        1000           0
PYRt2r              	       -1000          -0        1000           0
RPE                 	       -1000           0        1000           0
RPI                 	       -1000           0        1000           0
SUCCt2_2            	           0           0        1000        0.75
SUCCt3              	           0           0        1000           0
SUCDi               	           0           2        1000           0
SUCOAS              	       -1000          -2        1000           0
TALA                	       -1000           0        1000           0
THD2                	           0           0        1000         0.5
TKT1                	       -1000           0        1000           0
TKT2                	       -1000           0        1000           0
TPI                 	       -1000           1        1000           0</div></div></div></div><div class="inlineWrapper"><div  class = 'S12'><span style="white-space: normal"><span >ind=strcmp(model.rxns,</span><span style="color: rgb(170, 4, 249);">'EX_glc(e)'</span><span >);</span></span></div></div><div class="inlineWrapper outputs"><div  class = 'S10'><span style="white-space: normal"><span >FBAsolution_maxATP.rcost(ind)</span></span></div><div  class = 'S11'><div class='variableElement' style='font-family: Menlo, Monaco, Consolas, "Courier New", monospace; font-size: 12px; '>ans = 17.5000</div></div></div></div><div  class = 'S13'><span>Display the change in the flux vector:</span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S3'><span style="white-space: normal"><span >dv = FBAsolution_maxATP_moreGlc.v-FBAsolution_maxATP.v;</span></span></div></div><div class="inlineWrapper"><div  class = 'S4'><span style="white-space: normal"><span >dv(abs(dv)&lt;1e-4)=0;</span></span></div></div><div class="inlineWrapper outputs"><div  class = 'S10'><span style="white-space: normal"><span >printFluxVector(model, dv, 1)</span></span></div><div  class = 'S11'><div class="inlineElement eoOutputWrapper embeddedOutputsTextElement" uid="5FA2446E" data-testid="output_9" data-width="428" data-height="451" data-hashorizontaloverflow="false" style="width: 458px; max-height: 462px; white-space: normal; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"><div class="textElement" style="white-space: normal; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;">ACONTa              	           2
ACONTb              	           2
AKGDH               	           2
ATPM                	        17.5
ATPS4r              	        13.5
CO2t                	          -6
CS                  	           2
CYTBD               	          12
ENO                 	           2
EX_co2(e)           	           6
EX_glc(e)           	          -1
EX_h2o(e)           	           6
EX_o2(e)            	          -6
FBA                 	           1
FUM                 	           2
GAPD                	           2
GLCpts              	           1
H2Ot                	          -6
ICDHyr              	           2
MDH                 	           2
NADH16              	          10
NADTRHD             	           2
O2t                 	           6
PDH                 	           2
PFK                 	           1
PGI                 	           1
PGK                 	          -2
PGM                 	          -2
PYK                 	           1
SUCDi               	           2
SUCOAS              	          -2
TPI                 	           1</div></div></div></div></div><h4  class = 'S9'><span>What is your biochemical interpretation of this?</span></h4><h4  class = 'S9'><span>Hint: use drawFlux with a perturbed optimal reaction rate vector</span></h4><div  class = 'S2'><span>The flux map for the perturbation to optimal ATP production is shown below.  Note the reactions whose rates are substantially increasing, starting from glucose.</span></div><div class="CodeBlock"><div class="inlineWrapper outputs"><div  class = 'S14'><span style="white-space: normal"><span >drawFlux(map, model, dv, options);</span></span></div><div  class = 'S11'><div class="inlineElement eoOutputWrapper embeddedOutputsFigure" uid="5285C877" data-testid="output_10" style="width: 458px;"><div class="figureElement"><img class="figureImage figureContainingNode" src="" style="width: 560px;"></div></div></div></div></div><h4  class = 'S9'><span></span></h4><h4  class = 'S9'><span>Perturb the model in such a way as to increase the optimal rate of ATP hydrolysis ('ATPM') by exactly 17.5 units. How does this compare with the theoretical prediction?</span></h4><h4  class = 'S9'><span>Hint: change model.lb</span></h4><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S3'><span style="white-space: normal"><span >model = modelOri;</span></span></div></div><div class="inlineWrapper"><div  class = 'S4'><span style="white-space: normal"><span >model = changeRxnBounds(model,</span><span style="color: rgb(170, 4, 249);">'EX_glc(e)'</span><span >,-2,</span><span style="color: rgb(170, 4, 249);">'l'</span><span >); </span><span style="color: rgb(2, 128, 9);">%note the change in the lower bound from -1 to -2</span></span></div></div><div class="inlineWrapper"><div  class = 'S4'><span style="white-space: normal"><span >model = changeRxnBounds(model,</span><span style="color: rgb(170, 4, 249);">'EX_o2(e)'</span><span >,-1000,</span><span style="color: rgb(170, 4, 249);">'l'</span><span >);</span></span></div></div><div class="inlineWrapper"><div  class = 'S4'><span style="white-space: normal"><span >model = changeRxnBounds(model,</span><span style="color: rgb(170, 4, 249);">'ATPM'</span><span >,0,</span><span style="color: rgb(170, 4, 249);">'l'</span><span >);</span></span></div></div><div class="inlineWrapper"><div  class = 'S4'><span style="white-space: normal"><span >model = changeObjective(model,</span><span style="color: rgb(170, 4, 249);">'ATPM'</span><span >);</span></span></div></div><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: normal"><span >FBAsolution_maxATP_moreGlc = optimizeCbModel(model,</span><span style="color: rgb(170, 4, 249);">'max'</span><span >);</span></span></div></div></div><div  class = 'S13'><span>By changing the lower bound on glucose exhange from -1 to -2, we see that the value of the objective increases by 17.5, which is equal to the reduced cost of glucose obtained from FBAsolution_maxATP.rcost:</span></div><div class="CodeBlock"><div class="inlineWrapper outputs"><div  class = 'S14'><span style="white-space: normal"><span >FBAsolution_maxATP_moreGlc.f - FBAsolution_maxATP.f</span></span></div><div  class = 'S11'><div class='variableElement' style='font-family: Menlo, Monaco, Consolas, "Courier New", monospace; font-size: 12px; '>ans = 17.5000</div></div></div></div><h2  class = 'S1'><span>TROUBLESHOOTING</span></h2><div  class = 'S2'><span>Note that, if an optimization problem is reformulated from a maximisation to a minimisation problem, then the signs of each of the dual variables is reversed.</span></div><h2  class = 'S1'><span>TIMING</span></h2><div  class = 'S2'><span style=' font-style: italic;'>1 hr.</span></div><h2  class = 'S1'><span>ANTICIPATED RESULTS</span></h2><div  class = 'S2'><span>Understanding of how an optimal objective will change in response to changing the input data.</span></div><h2  class = 'S1' id = 'H_878295C9' ><span style=' font-style: italic;'>Acknowledgments</span></h2><div  class = 'S2'><span>Part of this tutorial was originally written by Jeff Orth and Ines Thiele for the publication "What is flux balance analysis?"</span></div><h2  class = 'S1'><span>REFERENCES</span></h2><div  class = 'S2'><span>1.    Orth. J., Thiele, I., Palsson, B.O., What is flux balance analysis? Nat Biotechnol. Mar; 28(3): 245–248 (2010).</span></div><div  class = 'S2'><span>2. Laurent Heirendt &amp; Sylvain Arreckx, Thomas Pfau, Sebastian N. Mendoza, Anne Richelle, Almut Heinken, Hulda S. Haraldsdottir, Jacek Wachowiak, Sarah M. Keating, Vanja Vlasov, Stefania Magnusdottir, Chiam Yu Ng, German Preciat, Alise Zagare, Siu H.J. Chan, Maike K. Aurich, Catherine M. Clancy, Jennifer Modamio, John T. Sauls, Alberto Noronha, Aarash Bordbar, Benjamin Cousins, Diana C. El Assal, Luis V. Valcarcel, Inigo Apaolaza, Susan Ghaderi, Masoud Ahookhosh, Marouen Ben Guebila, Andrejs Kostromins, Nicolas Sompairac, Hoai M. Le, Ding Ma, Yuekai Sun, Lin Wang, James T. Yurkovich, Miguel A.P. Oliveira, Phan T. Vuong, Lemmer P. El Assal, Inna Kuperstein, Andrei Zinovyev, H. Scott Hinton, William A. Bryant, Francisco J. Aragon Artacho, Francisco J. Planes, Egils Stalidzans, Alejandro Maass, Santosh Vempala, Michael Hucka, Michael A. Saunders, Costas D. Maranas, Nathan E. Lewis, Thomas Sauter, Bernhard Ø. Palsson, Ines Thiele, Ronan M.T. Fleming,</span><span> </span><span style=' font-weight: bold;'>Creation and analysis of biochemical constraint-based models: the COBRA Toolbox v3.0</span><span>, Nature Protocols, volume 14, pages 639–702, 2019</span><span> </span><a href = "https://doi.org/10.1038/s41596-018-0098-2"><span>doi.org/10.1038/s41596-018-0098-2</span></a><span>.</span></div><div  class = 'S2'></div>
<br>
<!-- 
##### SOURCE BEGIN #####
%% Sensitivity of a flux balance analysis solution with respect to input data
%% Author(s): Ronan M.T. Fleming, Leiden University
%% Reviewer(s): 
%% INTRODUCTION
% Consider an FBA problem
% 
% $$\begin{array}{ll}\textrm{max} & c^{T}v\\\text{s.t.} & Sv=b\\ & l\leq v\leq 
% u\end{array}\end{equation}$$
% 
% The local sensitivity of the optimal objective value $\mathcal{L}^{\star} 
% = c^{T}v^{\star}$  with respect to a changes in the input data $\left\lbrace 
% b,l,u\right\rbrace$ is given by
% 
% $$\frac{\partial\mathcal{L}^{\star}}{\partial b} = y^{\star} \\\frac{\partial\mathcal{\mathcal{L}^{\star}}}{\partial 
% l} = -w_{l}^{\star} \\\frac{\partial\mathcal{\mathcal{L}^{\star}}}{\partial 
% u} = w_{u}^{\star}$$
% 
% where $y^{\star}$ is a vector of shadow prices and $w = w_{l}-w_{u}$ is a 
% vector of reduced costs. That is, a shadow price is the partial derivative of 
% the optimal value of the objective function with respect to $b_i$.  It indicates 
% how much net production, or net consumption, of each metabolite increases (positive), 
% or decreases (negative), the optimal value of the objective. The reduced costs, 
% $-w_l$ and $w_{u\;}$are the partial derivative of the optimal value of the objective 
% function with respect to the lower and upper bounds on a reaction, respectively. 
% They indicate how much relaxation, or tightening, of each bound increases, or 
% decreases, the optimal objective, respectively. In the COBRA Toolbox, shadow 
% prices and reduced costs are calculated by |optimizeCbModel|. When using the 
% function
% 
% |FBAsolution = optimizeCbModel(model,'max');|
% 
% the shadow prices and reduced costs are given by |FBAsolution.y| and |FBAsolution.w|, 
% respectively. 
% 
% For a more complete theoretical description, see: cobratoolbox/tutorials/intro_sensitivityAnalysis.pdf
%% MATERIALS - EQUIPMENT SETUP
% Please ensure that all the required dependencies (e.g. , |git| and |curl|) 
% of The COBRA Toolbox have been properly installed by following the installation 
% guide <https://opencobra.github.io/cobratoolbox/stable/installation.html here>. 
% Please ensure that the COBRA Toolbox has been initialised (tutorial_initialize.mlx) 
% and verify that the pre-packaged LP and QP solvers are functional (tutorial_verify.mlx).
%% PROCEDURE
%% Load E. coli core model
% The most direct way to load a model into The COBRA Toolbox is to use the |readCbModel| 
% function. For example, to load a model from a MAT-file, you can simply use the 
% filename (with or without file extension). 

fileName = 'ecoli_core_model.mat';
if ~exist('modelOri','var')
modelOri = readCbModel(fileName);
end
%backward compatibility with primer requires relaxation of upper bound on
%ATPM
modelOri = changeRxnBounds(modelOri,'ATPM',1000,'u');
model = modelOri;
%setp the matlab e.coli metabolic map parameters
outputFormatOK = changeCbMapOutput('matlab');
map=readCbMap('ecoli_core_map');
options.zeroFluxWidth = 0.1;
options.rxnDirMultiplier = 10;
%% 
% 
% 
% The meaning of each field in a standard model is defined in the <https://github.com/opencobra/cobratoolbox/blob/master/docs/source/notes/COBRAModelFields.md 
% standard COBRA model field definition>.
% 
% In general, the following fields should always be present: 
%% 
% * *S*, the stoichiometric matrix
% * *mets*, the identifiers of the metabolites
% * *b*, Accumulation (positive) or depletion (negative) of the corresponding 
% metabolites. 0 Indicates no concentration change.
% * *csense*, indicator whether the b vector is a lower bound ('G'), upper bound 
% ('L'), or hard constraint 'E' for the metabolites.
% * *rxns*, the identifiers of the reactions
% * *lb*, the lower bounds of the reactions
% * *ub*, the upper bounds of the reactions
% * *c*, the linear objective
% * *genes*, the list of genes in your model 
% * *rules*, the Gene-protein-reaction rules in a computer readable format present 
% in your model.
% * *osenseStr*, the objective sense either |'max'| for maximisation or |'min'| 
% for minimisation
%% Sensitivity Analysis
% In the E. coli core model, when maximising ATP production, what is the shadow price of cytosolic protons? 
% Hint: |FBAsolution.y|

model = modelOri;
model = changeRxnBounds(model,'EX_glc(e)',-1,'l');
model = changeRxnBounds(model,'EX_o2(e)',-1000,'l');
model = changeRxnBounds(model,'ATPM',0,'l');
model = changeObjective(model,'ATPM');
printConstraints(model,-1000,1000)

FBAsolution_maxATP = optimizeCbModel(model,'max');
%% 
% Check the optimal value of the objective

FBAsolution_maxATP.f
%% 
% The shadow price of cytosolic protons (h[c]) is -0.25. 

ind=strcmp(model.mets,'h[c]');
FBAsolution_maxATP.y(ind)
printFluxVector(model, FBAsolution_maxATP.v, 1)
% What is your biochemical interpretation of this change in objective in the current context?
% Hint: printFluxVector, drawFlux
% This is a unique solution (see Example 3). 

dv = FBAsolution_maxATP_forceH.v-FBAsolution_maxATP.v;
dv(abs(dv)<1e-5)=0;
printFluxVector(model, dv, 1)
%% 
% The flux map for optimal ATP production is shown below.  

drawFlux(map, model, FBAsolution_maxATP.v, options);
%% 
% ATP production is constrained by cytoplasmic proton balancing.  Cytoplasmic 
% protons are produced by various metabolic reactions and also enter into the 
% cell, from the extracellular compartment, via the ATP synthase reaction (ATPS4r).  
% At steady-state, an equal number of protons must be pumped out of the cytoplasm 
% by the electron transport chain reactions or by excreting metabolites with symporters. 
% Setting model.b(i) = 4, where i corresponds to cytoplasmic protons, h[c], removes 
% 4 extra units of cytoplasmic protons from the system allowing 4 extra extracellular 
% protons to enter the system that then enter the cell via the ATP synthase reaction, 
% generating one extra unit of ATP. This increases the maximum rate of ATP synthesis 
% by one unit, thereby increasing the ATP yield from glucose by 1 mol ATP/mol 
% glucose. 
% Perturb the model in such a way as to increase the optimal rate of ATP hydrolysis ('ATPM') by exactly one unit. How does this compare with the theoretical prediction?
% Hint: change model.b
% Remove 4 units of cytoplasmic protons from the system, but changing model.b(i) 
% to 4, where i corresponds to the index for cytoplasmic protons, and calculate 
% the difference in the value of the optimal objective. The answer should be 1.

ind=strcmp(model.mets,'h[c]');
model.b(ind) = 4;
FBAsolution_maxATP_forceH = optimizeCbModel(model,'max');
FBAsolution_maxATP_forceH.f - FBAsolution_maxATP.f
% In the E. coli core model, when maximising ATP production, what is the reduced cost of glucose exchange? 
% Hint: FBAsolution.rcost

rcost = FBAsolution_maxATP.rcost;
rcost(abs(rcost)<1e-4)=0;
flux=FBAsolution_maxATP.v;
printFluxVector(model, [model.lb,flux,model.ub,rcost], 1)
ind=strcmp(model.rxns,'EX_glc(e)');
FBAsolution_maxATP.rcost(ind)
%% 
% Display the change in the flux vector:

dv = FBAsolution_maxATP_moreGlc.v-FBAsolution_maxATP.v;
dv(abs(dv)<1e-4)=0;
printFluxVector(model, dv, 1)
% What is your biochemical interpretation of this?
% Hint: use drawFlux with a perturbed optimal reaction rate vector
% The flux map for the perturbation to optimal ATP production is shown below.  
% Note the reactions whose rates are substantially increasing, starting from glucose.

drawFlux(map, model, dv, options);
% 
% Perturb the model in such a way as to increase the optimal rate of ATP hydrolysis ('ATPM') by exactly 17.5 units. How does this compare with the theoretical prediction?
% Hint: change model.lb

model = modelOri;
model = changeRxnBounds(model,'EX_glc(e)',-2,'l'); %note the change in the lower bound from -1 to -2
model = changeRxnBounds(model,'EX_o2(e)',-1000,'l');
model = changeRxnBounds(model,'ATPM',0,'l');
model = changeObjective(model,'ATPM');
FBAsolution_maxATP_moreGlc = optimizeCbModel(model,'max');
%% 
% By changing the lower bound on glucose exhange from -1 to -2, we see that 
% the value of the objective increases by 17.5, which is equal to the reduced 
% cost of glucose obtained from FBAsolution_maxATP.rcost:

FBAsolution_maxATP_moreGlc.f - FBAsolution_maxATP.f
%% TROUBLESHOOTING
% Note that, if an optimization problem is reformulated from a maximisation 
% to a minimisation problem, then the signs of each of the dual variables is reversed.
%% TIMING
% _1 hr._
%% ANTICIPATED RESULTS
% Understanding of how an optimal objective will change in response to changing 
% the input data.
%% _Acknowledgments_
% Part of this tutorial was originally written by Jeff Orth and Ines Thiele 
% for the publication "What is flux balance analysis?"
%% REFERENCES
% 1.    Orth. J., Thiele, I., Palsson, B.O., What is flux balance analysis? 
% Nat Biotechnol. Mar; 28(3): 245–248 (2010).
% 
% 2. Laurent Heirendt & Sylvain Arreckx, Thomas Pfau, Sebastian N. Mendoza, 
% Anne Richelle, Almut Heinken, Hulda S. Haraldsdottir, Jacek Wachowiak, Sarah 
% M. Keating, Vanja Vlasov, Stefania Magnusdottir, Chiam Yu Ng, German Preciat, 
% Alise Zagare, Siu H.J. Chan, Maike K. Aurich, Catherine M. Clancy, Jennifer 
% Modamio, John T. Sauls, Alberto Noronha, Aarash Bordbar, Benjamin Cousins, Diana 
% C. El Assal, Luis V. Valcarcel, Inigo Apaolaza, Susan Ghaderi, Masoud Ahookhosh, 
% Marouen Ben Guebila, Andrejs Kostromins, Nicolas Sompairac, Hoai M. Le, Ding 
% Ma, Yuekai Sun, Lin Wang, James T. Yurkovich, Miguel A.P. Oliveira, Phan T. 
% Vuong, Lemmer P. El Assal, Inna Kuperstein, Andrei Zinovyev, H. Scott Hinton, 
% William A. Bryant, Francisco J. Aragon Artacho, Francisco J. Planes, Egils Stalidzans, 
% Alejandro Maass, Santosh Vempala, Michael Hucka, Michael A. Saunders, Costas 
% D. Maranas, Nathan E. Lewis, Thomas Sauter, Bernhard Ø. Palsson, Ines Thiele, 
% Ronan M.T. Fleming, *Creation and analysis of biochemical constraint-based models: 
% the COBRA Toolbox v3.0*, Nature Protocols, volume 14, pages 639–702, 2019 <https://doi.org/10.1038/s41596-018-0098-2 
% doi.org/10.1038/s41596-018-0098-2>.
% 
%
##### SOURCE END #####
-->
</div></body></html>